###
#
# This replication script is for all analyses and figures associated with the manuscript:
# "Measuring Accountability in Interlocal Agreements between Indigenous and Local Governments"
#
# For additional details, please consult the README file included in the replication folder.
#
# Tyler Girard
# E: tgirard@purdue.edu
#
#
###

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Maintain Folder Structure Using "Here" Package ####

library(here)
here()

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Load Packages ####

library(rio)
library(tidyverse)
library(table1)
library(brms)
library(tidybayes)
library(ggpubr)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Prep BRMS ####

options(mc.cores = parallel::detectCores())

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Set Seed ####

set.seed(515)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Load Original Data ####

## Agreement Data

df <- rio::import("Muni-Indig Agreements (Replication File).xlsx") %>%
  mutate(id = 1:n()) %>%
  select(id, everything())

head(df)

ymat <- df %>%
  select(contains("ind_"))

ymat <- ymat*2

head(ymat)


## Load recode data

rc <- rio::import("Policy Area ReCode.xlsx")

rc <- rc %>%
  rename(PolicyArea = Original,
         PolicyArea_RC = New)

df <- df %>%
  left_join(rc)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Clean Data for 2PL-IRT Measurement Model ####

scaledat <- ymat %>%
  mutate(id = 1:n()) %>%
  pivot_longer(cols = -c("id"), 
               names_to = "item",
               values_to = "resp") %>%
  filter(item != "ind_WebSpace" &
           item != "ind_DollarValue") %>%
  mutate(resp = ifelse(resp > 1, 1, 0))


head(scaledat)

scaledat_wide <- scaledat %>%
  pivot_wider(names_from = item, values_from = resp) %>%
  select(-c(id)) %>%
  mutate(id = 1:n()) %>%
  select(c(id, everything()))

head(scaledat_wide)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Estimate Measurement Model ####

## Model Formula ##

formula_va_2pl <- bf(
  resp ~ exp(logalpha) * eta,
  eta ~ 1 + (1 |i| item) + (1 | id),
  logalpha ~ 1 + (1 |i| item),
  nl = TRUE
)

## Priors ##

prior_va_2pl <-
  prior("normal(0, 5)", class = "b", nlpar = "eta") +
  prior("normal(0, 5)", class = "b", nlpar = "logalpha") +
  prior("constant(1)", class = "sd", group = "id", nlpar = "eta") +
  prior("normal(0, 5)", class = "sd", group = "item", nlpar = "eta") +
  prior("normal(0, 5)", class = "sd", group = "item", nlpar = "logalpha") +
  prior("lkj(2)", class = "cor")

## Fit Model ##

start_m <- Sys.time()
modm <- brm(
  formula = formula_va_2pl,
  data = scaledat,
  family = brmsfamily("bernoulli", "logit"),
  prior = prior_va_2pl,
  chains = 4,
  seed = 2,
  cores = 4,
  control = list(adapt_delta = 0.99, max_treedepth = 15),
  file = "finalfits/modm",
  file_refit = "on_change"
)
end_m <- Sys.time()

end_m - start_m
# measurement model run time: ~17 mins


## Model Diagnostics ##

# Check 1: Divergences
# no divergences

# Check 2: Rhat
summary(modm)$random
summary(modm)$fixed
# all Rhats <= 1.01

# Check 3: Traceplots
pdf(file = "SI_Measurement Trace Plots.pdf")
mcmc_plot(modm, type = "trace")
dev.off()
# all trace plots suggest good mixing

## Extract Posterior Samples ##

postsamps_measure <- as_draws_df(modm) 

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Measurement Model Results and Figures ####

#### ~ ICC Plots ~ ####

post <- postsamps_measure %>% 
  select(b_eta_Intercept, b_logalpha_Intercept, starts_with("r_item")) %>% 
  mutate(iter = 1:n()) %>% 
  pivot_longer(starts_with("r_item")) %>% 
  mutate(item      = str_extract(name, "\\ind_+.*]$"),
         parameter = ifelse(str_detect(name, "eta"), "xi", "logalpha")) %>% 
  select(-name) %>% 
  pivot_wider(names_from = parameter, values_from = value)

head(post)

post$item <- gsub(",Intercept]","", post$item)
post$item <- gsub("ind_","", post$item)

head(post)

post <- post %>% 
  expand(nesting(iter, b_eta_Intercept, b_logalpha_Intercept, item, xi, logalpha),
         theta = seq(from = -6, to = 6, length.out = 100)) %>% 
  # note the difference in the equation
  mutate(p = inv_logit_scaled(exp(b_logalpha_Intercept + logalpha) * (b_eta_Intercept + theta + xi))) %>% 
  group_by(theta, item) %>% 
  summarise(p = quantile(p, c(0.025, 0.5, 0.975)), q = c(0.025, 0.5, 0.975)) %>%
  mutate(qlabs = fct_recode(as.character(q),
                            "lwr" = "0.025",
                            "med" = "0.5",
                            "upr" = "0.975"))

iccdat <- post %>% 
  select(-c(q)) %>%
  pivot_wider(names_from = qlabs, values_from = p) %>%
  mutate(itemfull = fct_recode(item,
                               "Annual Budget \nRequirement" = "AnnualBudgetRequirement",
                               "Annual Report \nRequirement" = "AnnualReportRequirement",
                               "Dispute Resolution" = "DisputeResolution",
                               "Prescribed Budget \nStandards" = "PrescribedBudgetStandards",
                               "Private Meeting" = "PrivateMeeting",
                               "Publicly Available" = "PubliclyAvailable",
                               "Public Meeting" = "PublicMeeting",
                               "Public Reporting \nRequirements" = "PublicReportingRequirements",
                               "Specific Partner \nRoles" = "SpecificPartnerRoles"))



#### ~ Latent Estimates ~ ####

lats <- postsamps_measure %>% 
  select(starts_with("r_id")) %>% 
  mutate(iter = 1:n()) %>% 
  pivot_longer(cols = starts_with("r_id"), names_to = "treaty", values_to = "xi") %>% 
  mutate(id = extract_numeric(treaty)) %>%
  group_by(id) %>% 
  summarise(xi = quantile(xi, c(0.025, 0.25, 0.5, 0.75, 0.975)), q = c(0.025, 0.25, 0.5, 0.75, 0.975),
            std.dev = sd(xi),
            avg = mean(xi)) %>%
  mutate(qlabs = fct_recode(as.character(q),
                            "lwr95" = "0.025",
                            "lwr50" = "0.25",
                            "med" = "0.5",
                            "upr50" = "0.75",
                            "upr95" = "0.975")) %>%
  select(-c(q)) %>%
  pivot_wider(names_from = qlabs, values_from = xi)

#### ~ Figures ~ ####

# ICCs
plot_icc_bw <- ggplot(iccdat, aes(x = theta, y = med, ymin = lwr, ymax = upr)) +
  geom_line(size = 1.2) +
  geom_ribbon(color = "black", alpha=0.2, linetype = "dashed") +
  facet_wrap(~itemfull) +
  labs(
    #title = "ICCs for the 2PL",
    #   subtitle = "Each curve is based on the posterior median and 95% credible interval.", 
    x = expression(theta~('Latent Estimate of Agreement Accountability')),
    y = expression(italic(p)(y==1))) +
  theme_minimal() +
  theme(
    legend.position = "none"
  )

plot_icc_bw

# Latent Estimates
plot_treaties <- ggplot(lats, aes(x = reorder(id, med), y = med)) +
  geom_linerange(aes(ymin = lwr95, ymax = upr95)) +
  geom_linerange(aes(ymin = lwr50, ymax = upr50), size = 1) +
  geom_point(shape = 20, color = "white", fill = "white") +
  theme_minimal() +
  theme(
    axis.text.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank()
  ) +
  labs(
    #title = "Latent Estimates for Each FN Treaty",
    #subtitle = "Each estimate includes the posterior median and 95% credible interval.",
    x = "Municipal-Indigenous Agreement",
    y = expression(theta~('Latent Estimate of Agreement Accountability'))
  ) +
  coord_flip()

plot_treaties

# Combined figure
plot_measurementresults_combined <- ggarrange(plot_icc_bw, plot_treaties,
                                              labels = c("A", "B"))


ggsave(plot_measurementresults_combined, filename = "Main_Measurement Model ICCs and Latents.png", scale = 1.55)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Prep Data for Predictive Models ####

df_pred <- df %>%
  mutate(provname = substr(Agreement, start = 1, stop = 2)) %>%
  select(id, provname, PolicyArea_RC) %>%
  mutate(PolicyArea_RC = as.factor(PolicyArea_RC),
         PolicyArea_RC = fct_relevel(PolicyArea_RC, "Land"))

dfmod <- lats %>%
  pivot_wider() %>%
  left_join(df_pred)


## Create 100 datasets for predictive model with measurement error

set.seed(515)

ndatasets <- 100

brmsdatamultiple <- list()

tmp_FN_est <- NULL

for(i in 1:ndatasets){
  
  for(j in 1:nrow(dfmod)){
    tmp_FN_est[j] = rnorm(1, dfmod$avg[j], dfmod$std.dev[j])
  }
  
  brmsdatamultiple[[i]] <- data.frame(
    avg = tmp_FN_est,
    PolicyArea_RC = dfmod$PolicyArea_RC
  )
  
}


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Descriptive Statistics Table for Predictive Model Data ####


descriptivetab <- dfmod %>%
  ungroup() %>%
  select(med, PolicyArea_RC)

table1::label(descriptivetab$med)<-"Accountability (Posterior Median)"
table1::label(descriptivetab$PolicyArea_RC)<-"Policy Area"

destab_out <- table1::table1(~med+PolicyArea_RC, 
               data=descriptivetab, 
               topclass="Rtable1-zebra",
               rowlabelhead="Variable")

destab_out

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Estimate Predictive Model *With* Measurement Error ####


start_p1 <- Sys.time()
modp_merr <- 
  brm_multiple(
    data = brmsdatamultiple, 
    family = gaussian,
    formula = avg ~ PolicyArea_RC,
    prior = c(prior(exponential(1), class = sigma),
              prior(normal(0, 5), class = b),
              prior(normal(0, 5), class = Intercept)),
    iter = 2500, 
    warmup = 1000, 
    cores = 4, 
    chains = 4,
    seed = 15,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "finalfits/modp_merr")
end_p1 <- Sys.time()

end_p1 - start_p1

# model run time: ~13 mins


## Model Diagnostics ##

# Check 1: Divergences
# no divergences

# Check 2: Rhat
round(modp_merr$rhats, 2)
# all values <= 1.01


## Re-estimate the same model, do not combine chains, in order to
## calculate Bayesian R2 for regression tables

## See: #https://github.com/paul-buerkner/brms/issues/997

modp_merr_r2 <- 
  brm_multiple(
    data = brmsdatamultiple, 
    family = gaussian,
    formula = avg ~ PolicyArea_RC,
    prior = c(prior(exponential(1), class = sigma),
              prior(normal(0, 5), class = b),
              prior(normal(0, 5), class = Intercept)),
    iter = 2500, 
    warmup = 1000, 
    cores = 4, 
    chains = 4,
    seed = 15,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    combine = FALSE # note this line
)

save(modp_merr_r2, file = "Bayes R2 Mod2.RData")
#load("Bayes R2 Mod2.RData")

#pool_R2 <- function(mlist, probs = c(0.025, 0.975), robust = FALSE, ...) {
#  r2_post <- sapply(mlist, bayes_R2, summary = FALSE, ..., simplify = TRUE)
#  posterior_summary(c(r2_post), probs = probs, robust = robust)
#}

r2_post <- sapply(modp_merr_r2, bayes_R2, summary = FALSE, simplify = TRUE)

posterior_summary(c(r2_post), probs = c(0.025, 0.975), robust = TRUE)
#Estimate     Est.Error       Q2.5              Q97.5
#0.0747       0.0279         0.0303             0.1390

posterior_summary(c(r2_post), probs = c(0.025, 0.975), robust = FALSE)
#Estimate     Est.Error       Q2.5              Q97.5
#0.0773       0.0280         0.0303             0.1390


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Estimate Predictive Model *Without* Measurement Error ####

start_p2 <- Sys.time()
modp_nomerr <- 
  brm(
    data = dfmod, 
    family = gaussian,
    formula = med ~ PolicyArea_RC,
    prior = c(prior(exponential(1), class = sigma),
              prior(normal(0, 5), class = b),
              prior(normal(0, 5), class = Intercept)),
    iter = 2500, 
    warmup = 1000, 
    cores = 4, 
    chains = 4,
    seed = 15,
    control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "finalfits/modp_nomerr"
)
end_p2 <- Sys.time()

end_p2 - start_p2
# measurement model run time: ~30 seconds


## Model Diagnostics ##

# Check 1: Divergences
# no divergences

# Check 2: Rhat
summary(modp_nomerr)
# all Rhats <= 1.01

# Check 3: Traceplots
pdf(file = "SI_Predictive (No Error) Trace Plots.pdf")
mcmc_plot(modp_nomerr, type = "trace")
dev.off()
# all trace plots suggest good mixing


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Predictive Model Results ####

summary(modp_merr, robust = TRUE)
summary(modp_nomerr, robust = TRUE)

bayes_R2(modp_nomerr, robust = TRUE)
# Estimate  Est.Error       Q2.5     Q97.5
# 0.114     0.0299         0.0635    0.173

summary(modp_nomerr, robust = TRUE, prob = 0.9)

modp_merr_posts <- posterior_samples(modp_merr)
modp_nomerr_posts <- posterior_samples(modp_nomerr)


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Contrast Plots (Policy Area) ####

#### ~ With Measurement Error ~ ####

orig_contrasts_policy_merr <- modp_merr_posts %>%
  select(-c(b_Intercept, sigma, lp__)) %>%
  gather() %>%
  group_by(key) %>%
  median_qi(.width = c(.95)) %>%
  rename(y=value,
         ymin=.lower,
         ymax=.upper,
         cat=key) %>%
  mutate(cat = fct_recode(cat,
                          "Land - Conservation" = "b_PolicyArea_RCConservation",
                          "Land - Cooperation" = "b_PolicyArea_RCCooperation",
                          "Land - Emergency Services" = "b_PolicyArea_RCEmergencyServices",
                          "Land - Health" = "b_PolicyArea_RCHealth",
                          "Land - Libraries" = "b_PolicyArea_RCLibraries",
                          "Land - Municipal Services" = "b_PolicyArea_RCMunicipalServices",
                          "Land - Recreation" = "b_PolicyArea_RCRecreation",
                          "Land - Social Services" = "b_PolicyArea_RCSocialServices",
                          "Land - Transportation" = "b_PolicyArea_RCTransportation",
                          "Land - Waste" = "b_PolicyArea_RCWaste",
                          "Land - Water" = "b_PolicyArea_RCWater"))

head(orig_contrasts_policy_merr)




df_combos <- t(combn(names(modp_merr_posts)[grep("Policy", names(modp_merr_posts))],2))

ls <- list()

for(i in 1:nrow(df_combos)){
  tmp <- modp_merr_posts %>%
    select(df_combos[i,])
  
  colnames(tmp)<-gsub("b_PolicyArea_RC","",colnames(tmp))
  
  tmp$dif <- tmp[,1] - tmp[,2]
  
  res <- median_qi(tmp$dif, .width = c(.95))
  
  res$cat <- paste0(colnames(tmp)[1], " - ", colnames(tmp)[2])
  
  ls[[i]] <- res
}

contrasts_policy_merr <- do.call("rbind", ls) %>%
  rbind(orig_contrasts_policy_merr) %>%
  mutate(rslt = ifelse(ymin <=0 & ymax >=0, "95% Credible Interval Contains 0", "95% Credible Interval Does Not Contain 0"))


#### ~ Without Measurement Error ~ ####

orig_contrasts_policy_nomerr <- modp_nomerr_posts %>%
  select(-c(b_Intercept, sigma, lp__)) %>%
  gather() %>%
  group_by(key) %>%
  median_qi(.width = c(.95)) %>%
  rename(y=value,
         ymin=.lower,
         ymax=.upper,
         cat=key) %>%
  mutate(cat = fct_recode(cat,
                          "Land - Conservation" = "b_PolicyArea_RCConservation",
                          "Land - Cooperation" = "b_PolicyArea_RCCooperation",
                          "Land - Emergency Services" = "b_PolicyArea_RCEmergencyServices",
                          "Land - Health" = "b_PolicyArea_RCHealth",
                          "Land - Libraries" = "b_PolicyArea_RCLibraries",
                          "Land - Municipal Services" = "b_PolicyArea_RCMunicipalServices",
                          "Land - Recreation" = "b_PolicyArea_RCRecreation",
                          "Land - Social Services" = "b_PolicyArea_RCSocialServices",
                          "Land - Transportation" = "b_PolicyArea_RCTransportation",
                          "Land - Waste" = "b_PolicyArea_RCWaste",
                          "Land - Water" = "b_PolicyArea_RCWater"))

head(orig_contrasts_policy_nomerr)




df_combos <- t(combn(names(modp_nomerr_posts)[grep("Policy", names(modp_nomerr_posts))],2))

ls <- list()

for(i in 1:nrow(df_combos)){
  tmp <- modp_nomerr_posts %>%
    select(df_combos[i,])
  
  colnames(tmp)<-gsub("b_PolicyArea_RC","",colnames(tmp))
  
  tmp$dif <- tmp[,1] - tmp[,2]
  
  res <- median_qi(tmp$dif, .width = c(.95))
  
  res$cat <- paste0(colnames(tmp)[1], " - ", colnames(tmp)[2])
  
  ls[[i]] <- res
}

contrasts_policy_nomerr <- do.call("rbind", ls) %>%
  rbind(orig_contrasts_policy_nomerr) %>%
  mutate(rslt = ifelse(ymin <=0 & ymax >=0, "95% Credible Interval Contains 0", "95% Credible Interval Does Not Contain 0"))



#### ~ Figures ~ ####

contrasts_policy_merr$mod <- "With Measurement Error"
contrasts_policy_merr$plotlab <- paste(contrasts_policy_merr$mod, contrasts_policy_merr$rslt)
contrasts_policy_nomerr$mod <- "Without Measurement Error"
contrasts_policy_nomerr$plotlab <- paste(contrasts_policy_nomerr$mod, contrasts_policy_nomerr$rslt)

contrasts_policy_combined <- rbind(contrasts_policy_merr, contrasts_policy_nomerr) %>%
  arrange(cat)

contrasts_policy_combined1 <- contrasts_policy_combined[1:66,]
contrasts_policy_combined2 <- contrasts_policy_combined[67:132,]

plot_contrasts_policy1 <- contrasts_policy_combined1 %>%
  ggplot(aes(y = y, x = cat)) +
  geom_hline(yintercept = 0, color = "indianred") +
  geom_point(aes(shape = plotlab, color = plotlab), position = position_dodge(width = 0.4)) +
  geom_linerange(aes(ymin = ymin, ymax = ymax, linetype = plotlab, color = plotlab), position = position_dodge(width = 0.4)) +
  #scale_color_manual(values = c("black", "dodgerblue")) +
  scale_colour_manual(name = "Results",
                      #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                      values = c("black", "black", "dodgerblue")) +   
  scale_shape_manual(name = "Results",
                     #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                     values = c(16,17,17)) +
  scale_linetype_manual(name = "Results",
                        #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                        values = c(1, 2, 2)) +
  theme_bw() +
  theme(
    legend.position = "right"
  ) +
  labs(
    x = "",
    y = "Estimate"
  ) +
  coord_flip()

plot_contrasts_policy1


plot_contrasts_policy2 <- contrasts_policy_combined2 %>%
  ggplot(aes(y = y, x = cat)) +
  geom_hline(yintercept = 0, color = "indianred") +
  geom_point(aes(shape = plotlab, color = plotlab), position = position_dodge(width = 0.4)) +
  geom_linerange(aes(ymin = ymin, ymax = ymax, linetype = plotlab, color = plotlab), position = position_dodge(width = 0.4)) +
  #scale_color_manual(values = c("black", "dodgerblue")) +
  scale_colour_manual(name = "Results",
                      #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                      values = c("black", "black", "dodgerblue")) +   
  scale_shape_manual(name = "Results",
                     #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                     values = c(16,17,17)) +
  scale_linetype_manual(name = "Results",
                        #labels = c("Control, Non-F", "Control, Flwr", "Exclosure, Non-F", "Exclosure, Flwr"),
                        values = c(1, 2, 2)) +
  theme_bw() +
  theme(
    legend.position = "right"
  ) +
  labs(
    x = "",
    y = "Estimate"
  ) +
  coord_flip()

plot_contrasts_policy2


# Save figures
ggsave(plot_contrasts_policy1, filename = "SI_Policy Contrasts 1.png", scale = 0.85)
ggsave(plot_contrasts_policy2, filename = "SI_Policy Contrasts 2.png", scale = 0.85)


#### ! Backup Save Files ! ####

save(modm, file = here("backups", "Measurement Model.RData"))
save(modp_merr, file = here("backups", "Predictive Model with M Error.RData"))
save(modp_nomerr, file = here("backups", "Predictive Model without M Error.RData"))
save(modp_merr_r2, file = here("backups", "Predictive Model with M Error for Bayes R2.RData"))

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#                                                     End of Main Script                                                              #
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#### Extra - Generate Latex Code for Linear Equation in SI File ####

library(equatiomatic)

polarea <- c("Land", "Conservation", "Cooperation", "Emergency Services", "Health", 
             "Libraries", "Municipal Services", "Recreation", "Social Services", "Transportation", "Waste", "Water")

equ_df <- data.frame(
  `Treaty Accountability` = rnorm(1000),
  `Policy Area` = as.factor(sample(polarea, 1000, replace = TRUE))
)

equ_df <- equ_df %>%
  mutate(Policy.Area = relevel(Policy.Area, "Land"))

equ_mod <- lm(Treaty.Accountability ~ Policy.Area, equ_df)

equatiomatic::extract_eq(equ_mod, wrap = TRUE, terms_per_line = 3)

